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ABSTRACT 

We present N-body simulation calculations of the dependence of the power spec- 
trum of non-linear cosmological mass density fluctuations on the equation of state of 
the dark energy, w — p/p. At fixed linear theory power, increasing w leads to an in- 
crease in non-linear power, with the effect increasing with k. By k = 10 /iMpc -1 , 
a model with w = —0.75 has ~ 12% more power than a standard cosmologi- 
cal constant model (to = —1), while a model with w — —0.5 has ~ 33% extra 
power (at z = 0). The size of the effect increases with increasing dark energy frac- 
tion, and to a lesser extent increasing power spectrum normalization, but is insen- 
sitive to the power spectrum shape (the numbers above are for Cl m — 0.281 and 
as = 0.897). A code quantifying the non-linear effect of varying to, as a function 
of k, z, and other cosmological parameters, which should be accurate to a few per- 
cent for k<10 ft,Mpc _1 for models that fit the current observations, is available at 
http://www.cita.utoronto.ca/~pmcdonal/code.html This paper also serves as an ex- 
ample of a detailed exploration of the numerical convergence properties of ratios of 
power spectra for different models, which can be useful because some kinds of numeri- 
cal error cancel in a ratio. When precision calculations based on numerical simulations 
arc needed for many different models, efficiency may be gained by breaking the prob- 
lem into a calculation of the absolute prediction at a central point, and calculations 
of the relative change in the prediction with model parameters. 



1 INTRODUCTION 



Dark energy is the most important focus of study in cos- 
mology today. Its basic existence is well established, both 
by observations of Type la supernovae jKnop et al.ll2003t 
iRiess et al.ll2004lh and by combinations of other observables 
(Seli ak et al.ll2005r> . The focus now is on probing the prop- 
erties of the dark energy (or whatever new physics causes 
the Universe to look like it contains dark energy). A simple 
first parameterization of its properties is to specify the ra- 
tio of pressure to density (equation of state) w = p/p, with 
w = — 1 for a cosmological constant. 

The large amount of observational effort focused on 
measuring w to high precision must be matched by suffi- 
ciently accurate predictions of the dependence of observ- 
ables on w. Theory calculations will need to be done more 
carefully than they have been in the past. In particular, 
probes of w based on cosmological structure (e.g., weak 
lensing, galaxy clusters, and even large-scale galaxy cluster- 
ing where the l inear power assumption is no longer exclu- 
sively relied on JPregmark et alJl20oj : lAbazaiian et alJl2004i 



lEisenstein et al.H2005lv rwill generally require non-linear nu 
merical simulations as the fundamental method of calculat 



ing theory predictions. This paper tackles a small part of the 
problem: computing the dependence of the non-linear mass 
power spectrum, Pt^i J (k,z), on w. 

The most direct use of calculations of the non- 
linear power spectrum is for weak gravitational lens- 
ing (c o smic shear) studies (e.g. , | 3enabed fc Van Waerbekel 
<2003y Simpson fc Bridle] (120041) : Huterer fc Takadal i2005ft : 
iJarvis et alJ J2005I) : iKnox et alJ <2005ft h Direct use of the 
mass power for weak lensing (as opposed to ray tracing) is 
of cours e only as good as th e Born and Limber approxi- 
mations. IVale fc White! (j2003) found the Born approxima- 
ti on agreed well with ra y tracing, although t his is less clear 
inlWhite fc Vald J2004f) . In a pilot project, IWhite fc Valel 
(2004) performed a small grid of simulations of weak lens- 
ing with full ray tracing, including w — —0.8, for parameter 
combinations designed to leave the CMB fluctuations in- 
variant. When performing a full grid of models for precision 
parameter fitting, this CMB-guided method is probably the 
best way to go. Here we concentrate on isolating the effect 
of w by running simulations where only w and possibly one 
other parameter is varied at a time (although the fitting 
formula we present does include joint variations of w, as, 
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and f2 m ). Even if ray tracing is ultimately required for high 
precision, it should still be useful to have an accurate mass 
power spectrum calculation, so we leave ray tracing for fu- 
ture work. 

The fitting formulas most commonly used to predict 
the non-linear power spectrum given a set of cosmologi- 
cal p arameters were not calibrated with w ^ — 1 simula- 
tions dPeacock fc Doddslll99el : ISmith et aljEooS) . although 
lBenabed^^BeTn £tfdeauT 2001 ) present an untested prescrip- 
tion for using IPeacock^^Doddsl dl996ll . Often these for- 
mulas are used for m / -1 by making the untested as- 
sumption that models with equal linear theory power and 
other parameters (most importantly Q m ) at the redshift of 
int erest will have equal non-linear power, independent of 
w. iMa et all (I1999T) did simulate the mass power spectrum 

3 , and combined these with 



for w = -2/3, -1/2, and -1 

the ACDM simulations in iMal Jl99ST) to produce a fitting 
formula. There is no obvious reason not to trust this for- 
mula to the 10% level of accuracy advertised; however, the 
grid of simulations used to calibrate it was sparse (e.g., the 
w > — 1 simulations all had Q. m > 0.4), and numerical con- 
vergence was not rigorously demonstrated. We will find that 
this formula does not work well in the region of par ameter 
space where we find ourselves. iKlypin et al.l i2003T) simu- 
lated models with different values of w, but did not present 
the power spectrum in detail. The one case they did show, a 
iRatra fc Peebles! il98ST> model with time varying w ~ —0.5, 
appears to be roughly consistent with our results. 

We do not aim in this paper to replace the existing 
fitting formulas for the non-linear power - only to provide 
an accurate correction for w ^ — 1. We generally consider 
flat models with parameters erg (the rms linear mass den- 
sity fluctuations in 8/i _1 Mpc radius spheres at z = 0), 
Q m , Qb (the mass and baryon densities relative to the crit- 
ical density), h (the Hubble parameter), n (the logarith- 
mic slope of the primordial power spectrum), and w. The 
dark energy equation of state parameter w does not affect 
the shape of the transfer function on non-linear scales (e.g., 
< 0.5% difference at fc > 0.007 ft Mpc - between the trans- 
fer function for w = -1 an d w = -1/2, from CMBFAST 
JSeliak fc Zaldarriaeal Il99rj) . for a typical model normal- 
ized at high-fc), so the z = linear theory power is not 
changed at all by changing w, when our other parameters 
are fixed. There can be, however, a change in the non-linear 
power with w, because the growth history changes - increas- 
ing w means the Universe had relatively more linear theory 
power at earlier times. This effect is not included in an esti- 
mate of the non-lin ear power made using something like the 
ISmith et all J2003I) fitting function. 

We will frequently refer to the change in observables 
with w at fixed z = values of other parameters as "the 
effect of u>." We are not implying that there is anything 
uniquely correct about this. From some points of view it 
would make more sense to fix the other parameters at early 
times, before the dark energy has become significant. We 
focus on the former definition of the effect of w simply 
because i t is generally not included in weak lensing cal- 
culations dHuterer fc Takadall2005t jSimpson fc Bridldl2004 
iJarvis et alJl2005l iKnox et al.ll2005l) ." and can not be repro- 
duced in any obvious way using the existing codes. While it 
is mostly a matter of arbitrary definition, there is a real sense 



in which our choice is "the" non-linear effect of w at z — 0: 
our effect goes to zero in the limit of small perturbations. 

The range of fc and z we consider, and the accuracy 
goal, are guided by the weak lensing application (fc = 2tt/\, 
where A is the wavelength of a Fourier mode). We limit 
ourselves to the range < z < 1.5, because this range 
is most directly sensitive to the presence of dark energy 
and most relevant to galaxy weak lensing surveys, and be- 
cause the power at increasing redshift should probably be 
simulated using decreasing box size, since limited particle 
density becomes an increasingly serious problem (the true 
small-scale power is smaller relative to the spurious par- 
ticle discreteness-related power) while limited box size be- 
comes less problematic (s maller scale modes are still linear) . 
iHuterer fc Takadal i2005T) investigated the requirements on 
the PNh(k, z) calculation for future large weak lensing sur- 
veys, finding that 1-2% accuracy should be sufficient, or 
0.5% in the worst possible case. Th ey found that the rel- 
evant scales are 0.1<fc<10 h Mpc" 1 . Izhan fc Knoxl (120041) 
studied the effect of hot baryons on the weak-lensing shear 
power spectrum in halo models. They found an effect of 
roughly 5-10% on the power at fc = 10 h Mpc -1 (reading 
from their Figure 1), so trying to achieve better than a few 
percent precision at this fc using simulations without gas 
dynamics is probably pointless (it is always good to aim 
for er r ors so mewhat smaller than the other known sources) . 
IZhanl (120041) found roughly similar results i n SPH simula- 
tions (reading from Figure 5.7), and I White! <l2004h found a 
comparable scale for the effect of baryonic cooling. We will 
not actually achieve 1-2% level accuracy in the dependence 
of power on w, but we think we show a clear path to it. The 
accuracy we do achieve is substantially better than anything 
else available, so it should be useful until a larger project can 
improve it. 

Throughout the paper, we use the trick of canceling nu- 
merical errors by taking ratios of power spectra for different 
models, rather than looking at the absolute power in each 
model directly. This should in no way be considered a swin- 
dle or an added approximation, ft is completely reasonable 
and expected that many types of error are insensitive to the 
model, so that ratios (or differences) between the models can 
genuinely be computed more accurately than either model 
individually. One example of this is the effective smooth- 
ing involved in mapping particles to a grid so that you can 
FFT the periodic density field for the purpose of measuring 
the power. This suppresses the power substantially, but by 
precisely the same factor in all models. This factor cancels 
exactly when we take a ratio, up to some high fc where the 
power may be corrupted by aliasing, or the suppression is 
simply too large to invert accurately. The key to believing 
these results is that the convergence of the ratios of power 
with numerical parameters of the simulations must be tested 
carefully in the same way a direct measurement of the power 
would be. For example, taking ratios would not cancel an 
additive Poisson noise term due to the limited number of 
particles, but this would become completely obvious in a 
convergence test where the number of particles is varied. 

We use a particle-multi-mesh (PMM) code, based on an 
improved particle-mesh (PM) algorithm, for our grid of N- 
body simulations. In principle, PM codes can achieve high 
spatial resolution but at a great cost in memory and to a 
lesser extent in work. In practice, they are normally lim- 
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ited to a mean interparticle spacing-to-mesh cell spacing ra- 
tio of 2:1, where the storage requirements for particles and 
grids are approximately balanced. PMM utilizes a domain- 
decomposed, FF T-based gravity solver dTrac fc Penll2004 
iMerz et al.ll2005F) to achieve higher spatial resolution while 
maintaining memory costs. We will find that a ratio of 4:1 
between the particle grid spacing and mesh grid spacing is 
roughly optimal, in the sense that at a given k (the most rel- 
evant k where the errors are a few percent) the error from 
limited particle density and limited force resolution are sim- 
ilar, for ratios of pow er spectra from s i mulat ions with differ- 
ent w. We note that iHeitmann et alJ ((2004) compared sev- 
eral N-body codes (not including ours), looking at absolute 
power, and found good general agreement. 

It can be useful to consider numerical errors rather care- 
fully, focusing on the initial small breakdown of accuracy, be- 
cause, for example, a factor of two change in resolved scale 
or required box size can easily change the require computer 
time or memory by a factor of at least eight. At least five 
things must be investigated in every simulation program: 
Sensitivity to starting redshift, box size, mass resolution 
(i.e., number of particles), force resolution, and time step 
size. We also test our method of power spectrum calcula- 
tion. Even though not everything we find will be completely 
generalizable to other types of simulations, we think it is 
useful to lay out an attempt to push them all to a percent 
level of control. 

Throughout this paper we use the philosophy of sup- 
pressing our urge to be general (e.g., increase the redshift 
range and parameter space coverage and study weak lens- 
ing directly) in favor of patiently investigating a restricted 
problem. I.e., rather than trying to entirely solve the prob- 
lem of calculating the power spectrum over all parameter 
space, which at this point would require approximations and 
cutting corners, we select a small subproblem and explore 
it carefully, in hopes of learning better how to do the full 
problem accu rately and effi c iently . 

Recently iKuhlen et alJ !|2()0~> i performed a set of sim- 
ulations of models with different w values, although they 
focused on dark matter halo density profiles. We mention 
them because they also represent a good source for some 
pedagogical figures, e.g., of the linear transfer function and 
growth factor as a function of w. 

The rest of the paper is laid out as follows: First, in 
§2, we qualitatively demonstrate the effect of changing the 
equation of state of the dark energy on the non-linear mass 
power spectrum. Then, in §3, we describe detailed tests of 
our simulations (this section is not intended for the casual 
reader). Finally, in §4, we describe the code we provide to 
quantify the results. 

We often describe simulations using the notation 
(L, P, M), where L is the box size in comoving h~ 1 Mpc, 
P 3 is the number of particles, and M 3 is the number of mesh 
cells. 



2 THE EFFECT OF CHANGING W 

Figure^shows the basic effect of w on the non-linear power, 
at z = 0. To minimize the statistical fluctuations, we have 
used the same initial conditions for simulations with differ- 
ent w, except that we adjust the initial amplitude of the 
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Figure 1. Fractional effect of w on the non-linear mass power 
spectrum, at fixed linear theory power. The thick lines of a given 
color/type show, from top to bottom, w = —0.5, —0.75, —1.5, all 
relative to to = —1. The black, solid line is for Q m = 0.281, with 
red, dashed (green, dotted) showing f2 m = 0.211 (0.351). Thin 
black lines show rms statistical error bands (in Figures IH4I the 
statistical errors are all similar to these examples). 



perturbations by the factor necessary to produce an iden- 
tical linear theory density field at z = (and as a result 
identical ag). All of the figures in this section show results 
from (110,192,768) simulations. In the next section we inves- 
tigate the accuracy of the results in detail, finding statistical 
and numerical convergence to a few percent. Inevitably, the 
effect of w increases with decreasing Sl m , i.e., increasing dark 
energy fraction, as shown in the figure. To be clear: in our 
many figures like this, we are plotting the ratio of power 
in a model with aj / -1 to a model with w = — 1, with 
all the other parameters the same at z = in both. For 
example, when we show the result for a different value of 
fi m , fim has been changed for both the model in the nu- 
merator and the model in the denominator For notational 
compactness, we will sometimes refer to this type of ratio as 
f w {k) = P w {k)/P w =-x{k). 

To put the 12% (33%) effect of w = -0.75 (w = -0.5) 
at k = 10 ft Mpc -1 in perspective, we note that the ra- 
tio of linear growth of power from early times to z — 
between the w = —1 and w — —0.75 (-0.5) models is 
D 2 (w)/D 2 (w = -1) = 0.83 (0.57). Note that weak lens- 
ing measur ements tend to be most sensitive to somewhat 
smaller k fauterer fc Takadal Eooil) . On the other hand, 
iHaean et alJ J2005T) recently showed that a modification of 
the power similar in form to ours has a substantial effect 
on the convergence power over a wide range in I. Consid- 
ering the non-trivial k a nd z dependence, a full p arameter 
forecast calculation fe.g.. Isimpson fc Bridle! (120041) ') will be 
needed to see if this effect can change the projections for 
future weak lensing measurements of w significantly. 
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Maetal 



We h ave c ompared this result to the formula of 
1999), and the agreement is not good. The 
1999) formula predicts a much larger effect of 
f w =-o.5(k) passing 1.5 at k = 0.4/iMpc -1 , be- 
2. This disagreement is not 
did not use simulations 



w, with f w =-o.s(k) passing 
fore plateauing at a value of ■ 
too surprising since iMa et alJ (11999 
in which w was varied at fixed values of the other param- 
eters. Without a very comprehensive grid of simulations, 
the simple dependence of the non-linear power on the lin- 
ear power (and Q m ) could easily be mixed with the type 
of w dependence that we are studying here. We remind the 
reader that no rigorous anal ytic prediction fo r t he non-linear 
power spectrum exists. The lMa et al.l 1119991) or lSmith et alJ 
( 2003) formulas are not derived from first principles - they 
are physically motivated fitting functions that are used to 
interpolate/extrapolate simulation results. Even if their ba- 
sic motivation is perfectly valid, they contain completely 
free parameters that are only determinable through fits to 
simulations, so their predictions are only as good as the sim- 
ulations used to calibrate them. IMa et ahl il999t) used only 
five models, with (w, erg, ii m , h) = (-1, 1.29, 0.3, 0.75), (- 
1, 1.53, 0.5, 0.7), (-2/3, ?, 0.4, 0.65), (-1/2, ?, 0.4, 0.65), 
and (-1/3, ?, 0.45, 0.65). The value of <rg was not given for 
the w > — 1 simulatio ns, but they were COBE normalized 
feunn fe Whitd|l997ll . meaning they all have different erg. 
One might hope that this was enough to calibrate the fitting 
formula everywhere, but this is not at all obvious. It would 
be very difficult using this set of simulations to disentangle 
the effects of the different parameters well enough to accu- 
rately extrapolate to a point in parameter space not bounded 
by the set. In particular, a parameter that has a relatively 
small direct effect will be most difficult to control. 

The normalization of the power spectrum, erg, also af- 
fects the result in a small, qualitatively reasonable way, with 
increasing w dependence for increasing as, as shown in Fig- 
ure [5] (the non- linear effect must disappear as as — > 0). 
Figure also shows that the changes in the shape of the 
power spectrum that arise from changes in h, and n do 
not change the w dependence significantly. 

To explore the origin of the w effect, in Figure|3]we show 
a similar calculation of the effect of changing Q m = 1 — Oa 
while holding the z = linear power fixed (along with 
all the other parameters, including w = — 1). Reducing 
Cl m to 0.192 (from 0.281) produces a model with the same 
power at 2 = 24 (our starting redshift) as the model with 
w — —0.75 and Q m = 0.281. We see that the non-linear 
power in these two models is quite similar. When we match 
the linear power at z = 1.9, using Q m — 0.213, the non- 
linear power is even more similar. It seems likely that Fig- 
ure and the effect of w in general, could be interpreted 
in a halo model by accounting for the difference in forma- 
tion redshift, and thus density profile, of the typical halos 
dominating the power on the scale of intere st (||^hlen et^tl] 
| 2005l: iHuffenberger fe Seliakl 120031: iBartelmann et alJ 12003 
bolae et al.ll2004l: iBartelmann et alj|200l l2002t) . 

The code of lSmith et al.l i2003l) should be able to repro- 
duce this fi m dependence. As we see in Figure |3 the agree- 
ment is excellent (we can not say for sure that our simula- 
tions are accurate enough to believe the ~ 5% disagreement 
at k = 10/iMpc -1 , although they probably ar e). Note that 
to ma ke this comparison we have modified the lSmith et alJ 
( 2003) code to accept an input linear theory power spectrum 
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Figure 2. Similar to Figure HI except thick red, dashed (green, 
dotted) lines show as = 0.800 (0.994). (our central model has 
trg = 0.897). Thin lines for variations U b = 0.0462 ± 0.0052, 
h = 0.710 ±0.066, and n = 0.980 ±0.065 are also plotted, to show 
that they are usually indistinguishable from the central model. 

in pla ce of its usual calc ulation base d onlBond fc Efs tathioul 
jl984lh The formula of E GUI); IMa et all l)l999h does 
poorly in this test, probably because its Q m dependence 
was essentially calibrated by two simulations with Q m — 0.3, 
as = 1.29, and fl m — 0.5, erg = 1.53, i.e., far from the com- 
bination of parameters we are testing. 

Finally, in Figure 0] we show the effect of w on the non- 
linear power at z — 1.5, at fixed values of our other param- 
eters. Note that in Figure 2] we are comparing models with 
different linear theory power at the observed redshift, and 
different £l m (z). The difference between Pnl(z = 1.5, w = 
—0.5) and Pnl{z = 1.5, w = —1) at fixed z — 1.5 linear 
theory power and Q m {z = 1.5) is small, so it is not very 
interesting to plot this, but we show one case, comparing 
w = —1.0 with our usual linear theory power spectrum at 
z = 0, but Qm = 0.150, to w = —0.5 with as increased 
to 1.004, and fi m = 0.411 (these two models have exactly 
matching linear power and Q m (z) at z = 1.5). 



3 NUMERICAL DETAILS 

In this section we investigate the dependence of our cal- 
culations on the numerical parameters of the simulations. 
Beyond testing the specific results we present, we hope to 
contribute to the collective wisdom of the research commu- 
nity about how to do precision cosmology based on numeri- 
cal simulations by exploring the idea of looking at the ratio 
of power in different models in simulations using identical 
random numbers for the initial conditions. 

For our grid of N-body simulations, we used the PMM 
code, an improved version of the PM algorithm. It is based 
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Figure 3. Comparison of the effect of changing w to the effect of 
changing f2 m , at fixed z = linear theory power. Black, solid line: 
our standard ratio of simulated Pnl(^) with w = —0.75 to w = 
— 1. Blue, long-dashed line: ratio of power with Q m = 0.192 to our 
standard Q m = 0.281 (both with w = —1). These two alternative 
models have identical linear theory power at both z = 24 and 
2 = 0. Red, short-dashed line: ratio of power with Q m = 0.213 
to standard (this case has linear power equal to the ui = —0.75 
case at z = 1.9). Green, dotted line: as blue, but based on Smith 
et al. (2003) (for both numerator and denominator). Magenta, 
dot-dashed: as blue but based on Ma et al. (1999). 



on a two-level mesh Poisson solver where the gravitational 
forces are separated into long-ran ge and short-range com- 
ponents, as descr ibed in detail in Trac fc Pen! 1^001) and 
iMerz et al.l <2005l) . The long -range force is computed on the 
root-level, global mesh, much like in a PM code. To achieve 
higher spatial resolution, the domain is decomposed into cu- 
bical regions and the short-range force is computed on a 
refinement-level, local mesh. In the current version, PMM 
can achieve a spatial resolution of 4 times better than a 
standard PM code at the same cost in memory. 

Simulations with P = 192 and Al — 768 fit conveniently 
into a single node on the 528-CPU Beowulf cluster at CITA, 
and run in less than a day. The importance when doing 
this kind of study of complete freedom to run an arbitrar- 
ily large number of exploratory simulations, with relatively 
quick turn-around, can not be underestimated, so we focus 
on this configuration. Future precision simulation grids will 
of course use larger simulations. Our initial guess was that 
(220, 192, 768) simulations would be most useful, so this sec- 
tion will focus first on the convergence properties of these, 
including comparisons of (110, 96, 384) to (110, 192, 384), 
to test the effect of finite particle density at the same force 
resolution as the (220, 192, 768) simulations, and compar- 
isons of (110, 96, 384) to (110, 96, 768) to test the effect of 
force resolution. Ultimately, we will conclude that we were 
overly worried about limited box size and insufficiently wor- 



Figure 4. Non-linear power for varying w at z = 1.5. Thick 
black (solid) lines are for models with identical linear power and 
f2 m at z = (from top to bottom w = —0.5, -0.75, and -1.5, 
all relative to ui = —1). Thin black lines show the rms statistical 
error bands. Note that much of the change here is accounted for by 
simple differences in linear growth. The red (dotted) line shows 
w = —0.5 relative to w = — 1 when the two models have been 
constructed to have identical linear power and Q, m (z) at z = 1.5. 



ried about limited resolution, so subject to the constraint 
P = 192, M = 768, somewhat smaller box size is optimal 
(the code we release is based on (110, 192, 768) simulations). 

Throughout this section our plots use a standard verti- 
cal axis range 0.94 — 1.06, to allow easier comparison of the 
size of different errors. 



3.1 Initial conditions 

Our transfer functions are computed using the " l ingers " 
code associated with grafic2-1.01 feertschingerl 1200 ll) . 
GRAFIC2 is then used to generate the initial conditions. We 
turn off the GRAFIC2 Hanning window, which isotropizes 
the small-scale structure at the expense of suppressing the 
small-scale power. 

To determine the linear growth factor (used to convert 
initial conditions generated for w = — 1 models into w / -1 
models) we numerically solve 



D" + l 



w(a) 



D' 3 X(a) D 



1 + X(a) 
finder fc .TenkinsiEoO^. with 



2 l + X(a) a 2 



: 



X(a) = 



- ^ n 



l - fi„ 



-3 T„ d In a' w (a' ) 



(i) 



(2) 



When we compare simulations with the same box size 
but different particle density, the initial conditions of the box 
with fewer particles are set by a sharp fc-space filter applied 
to the initial conditions of the higher resolution box. This is 
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Figure 5. Change in the effect of w, f w ——o.s(k) = 
P w =—0.5{k)/ P w =-i, with starting redshift in the simulations. 
Thick (thin) lines are power at z = (z = 1.5). The denominator 
is for starting z\ = 799, with black/solid showing the difference 
for zi = 24, blue/dotted z; = 49, green/short-dashed Z[ = 99, 
magenta/long-dashed z\ = 199, and red/dot-dashed z\ = 399. 
The two black lines of each type show different realizations of the 
initial conditions. 



equivalent to simply regenerating the initial conditions with 
the same random numbers for the large-scale modes (for our 
method of generating initial conditions), but not equivalent 
to re-binning the density and momentum fields in real space 
(the latter method introduces a suppression of high-fc power 
that increases the level of disagreement between the results 
for the two particle densities). 



3.2 Starting redshift 

It is important for precision calculations to test the affect 
of changing the starting redshift in the simulations. This 
can identif y, for example, tr ansients due to t he imperfec- 
tion of the IZel'Dovichl dl970l) approximation dScoccimarrol 
ll99Sf) . Figure |S| shows the change in the effect of w as we 
increase the starting redshift, Zi, from our standard Zi = 24. 
Note that it is not automatically the case that higher Zi is 
better, because numerical errors (most obviously suppres- 
sion of power by limited force resolution) have more time 
to accumulate in that case. The agreement is good but not 
perfect, with errors as large as 2% at z — and ~ 3% at 
z — 1.5. Subsequent to our runs for this paper, we found 
a small problem in the PMM force kernel which leads to 
this disagreement - it is not serious so we leave the more 
accurate calculation for the future. 



Figure 6. Change in the effect of w, f w ——o.s(k) = 
P w =— 0.5(h)/ Pw=—l, with number of particles. Thick (thin) lines 
are power at z = (z = 1.5), with L = 110 /i -1 Mpc (a) or 
L = 55h _1 Mpc (b) and force mesh M = 384. The numera- 
tor is simulations with (P = 96) 3 particles, the denominator 
is P = 192. Red (dashed) lines include subtraction of Poisson 
noise P no isc(k) = (L/P) 3 , black (solid) do not. For reference, 
the green (dotted) lines show 1 + P no ise(k) / Pme&sured{k) , where 
^measured (k) is the measured power for the w = — 1 model (af- 
ter deconvolution of the power spectrum measurement mesh, but 
without noise subtraction). The two green lines of each thickness 
represent the two particle densities. The vertical cyan (dotted) 
lines mark k = (96, 192)tt/110 h' 1 Mpc. 

3.3 Mass resolution 

Figure |SJ a) shows the effect of varying the number of parti- 
cles, for a fixed force resolution, by comparing the average 
of two (110, 96, 384) simulations to the average of two (110, 
192, 384) simulations (the pairs have different random initial 
conditions). Note that this is a test of the effect of adding 
high k power in the initial conditions in addition to the ef- 
fect of simply subdividing the mass. We have done the same 
comparison with a factor of two better force resolution and 
the results are similar. We see that for a (110, 96, 384) sim- 
ulation at z = 0, limited particle density becomes a >2% 
problem only at fc>8 h Mpc" 1 . At z = 1.5 it is a more seri- 
ous problem, surpassing 3% at k ~ 4 h Mpc -1 and quickly 
diverging (although note the expanded axis scale - the error 
is actually only 15% at k = 10 h Mpc -1 ). 

A potential solution to this problem is to subtract, after 
correcting for the mass assignment smoothing, Poisson shot- 
noise power -Pnoise ffc) = , where n = (P/L ) 3 is the mean 
parti cle density feaueh fc Efstathioul Il994l iB aueh et all 
1995 ^); however, as others hav e noted feaueh et al l Il995t 
ISmith et alll2003t ISirkol 120051) . the idea that the effect of 
finite particle density is to add this white noise component 
is only a guess, not something that can be assumed to hold, 
and in fact it is known not to hold at early times for a uni- 
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form grid start. Figure [5] calls into question the idea that 
subtracting white shot-noise i s ever useful for high precision 
calculations. (See Figure 11 of lSirkol l|2005f) for an enlighten- 
ing plot of the evolution of the particle discreteness power 
starti ng from a fixed grid - we produced a similar figure, 
but ISirkol j200^ 's is similar enough that it is not worth 
including ours in this paper.) Note that when we say that 
subtracting Poisson noise is not useful for high precision cal- 
culations, this does not mean it never leads to an improve- 
ment in accuracy - it sometimes does - the problem is that 
there is no clearly identifiable regime where the correction 
is both significant and accurate enough for high precision 
work. Th is should probably be regarded as a well-known 
fact (e.g.. lBaugh et al] lll995l) conclude that the discreteness 
corrections they discuss can not be applied consistently, and 
simply resort to using enough particles to make them neg- 
ligible), but i t is worth reiterat ing. It seems unlikely that 
a glass start llSmith et al.ll2003T) will lead to perfectly sta- 
ble, non-interacting, discreteness power either. Ultimately, 
direct tests of the convergence of observable statistics with 
increasing particle den sity for different methods of setting up 
the initial conditions iSmith et al]l2003| I Joyce et all 12004 
ISirkol200^1 . and possibly different methods for correcting for 
discreteness, are the only way to determine which method 
works best. 

Figure |SJb) shows the same comparison as Figure |SJb), 
reduced in scale by a factor of two, i.e., (55, 96, 384) is com- 
pared to (55, 192, 384), to estimate the accuracy of a (110, 
192, 768) simulation (there is some noise in this comparison 
so we have averaged two realizations of each size simulation). 
The results are much better, with accuracy better than 2% 
at z = 0, and better than 4% for z = 1.5 (better than 2% for 
fc<6 h Mpc -1 ). Note that the apparent helpfulness of Pois- 
son noise subtraction at z — 1.5 is probably coincidental, as 
it quickly becomes an over-correction at k > 10 /i Mpc -1 . 

Our default in this paper is to not subtract shot noise. 



Figure 7. Change in the effect of w, f w =—0.B{k) = 
P w —-0.5(k)/P w ——l, with force resolution. Black/solid 
(red/dotted) lines are power at z = (z = 1.5) in 
L = 110 h~ 1 Mpc simulations with (P = 96) 3 particles 
(thick lines) or P = 192 (thin lines). The denominator is for 
mesh M = 768, with numerator M = 384. The blue/long-dashed 
(green/short-dashed) lines show the same comparisons for 
L = 55 h~ 1 Mpc simulations (the effect appears to change sign, 
i.e., the curves are not mislabeled). For reference, magenta/dot- 
dashed lines show the ratio of power spectra [i.e., simply P(k), 
not fm(k)] for w = —1, comparing M = 384 to M = 768, for an 
L = 55 h. -1 Mpc box with P = 96 (thick) or P = 192 (thin) (the 
poorer agreement in each case is z = 0, better is z = 1.5). The 
vertical cyan/dotted lines show k = (96, 192)7r/110 h~ 1 Mpc. 

3.4 Force resolution 

Figure shows the effect of increasing the force mesh res- 
olution, for a fixed number of particles. Comparing L = 
110 /i -1 Mpc simulations with (M = 384) 3 force mesh cells 
to M = 768, we see that the error on the former is as large 
as 4% at z = and 7% at z = 1.5. The errors fall to < 3% 
for a similar comparison using L = 55 ft -1 Mpc boxes. In 
Figure |7] we break from our general policy of plotting only 
ratios of different models to show the lack of convergence of 
the absolute power spectrum, PM=3S4.(k) /Pm =768 (k). This 
shows the value of computing ratios - the better than 3% 
agreement in / ro =-o.s(fc) for the two meshes occurs in spite 
of a suppression of the absolute power by as much as 36%. 

3.5 Box size 

Insufficient box size can cause two problems: simple random 
fluctuations around the mean of a statistic due to limited 
volume, i.e., sample variance; and systematic errors in the 
mean of a statistic due to missing couplings to large-scale 
modes (or even small-scale modes missing due to limited 
fc-resolution) . The first of these can be eliminated by av- 
eraging over sufficient realizations of the initial conditions 
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Figure 8. Level of random fluctuations in the fractional effect 
of w on the non-linear power. Black, thin lines show fw=—o.s{k) 
from L = 220/i _1 Mpc simulations with nine different random 
seeds, each divided by the average of f w ——0.5(k) over all nine, 
at z = 0. The inner thick red/solid (green/dotted) lines show 
bin-by-bin the standard deviation of the nine around their mean 
(i.e., the error when using a single realization) at z = (z = 1.5). 
The outer thick lines are the same but for L = 110 /i -1 Mpc 
simulations. 
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Figure 9. Effect of box size on the fractional effect of w on the 
non-linear power. The black/solid (bluc/short-dashcd) line shows 
at z = (z = 1.5) the ratio of f w —— o.s(fc) computed by aver- 
aging over nine L = 110 /i — 1 Mpc simulations to the result from 
L = 220 /i — 1 Mpc simulations similarly averaged. The red/long- 
dashed (green/dotted) line shows at z = (z = 1.5) similar ratios 
for L = 220 ft -1 Mpc compared to L = 440 h -1 Mpc. Note that 
different sized simulations do not have matching grids in k, so 
these comparisons involve some interpolation. 



while the second can not (although various methods ha ve 
been proposed to i mprove the results, e.g., ISirkol l)200fi|) ;i. 
iBaela fc Ravi <l2005l) discuss requirements on simulation box 
size, but not for the power spectrum - note that the require- 
ments on numerical parameters will generally be different for 
different statistics (one advantage of the power spectrum is 
that it is not directly sensitive to structure on scales larger 
than the box). 

Figure [S] shows that a single L — 220 ft -1 Mpc sim- 
ulation is sufficient to compute the fractional difference in 
power between w = —0.5 and w = — 1 to better than 1% 
rms statistical error. A single L — 110 h~ Mpc simula- 
tion would achieve about 2% precision at z — and 3% 
at z — 1.5. In principle this Figure can be sensitive to the 
binning in k but in practice it is not because the errors in 
nearby bins are strongly correlated. 

Figure |5] addresses the systematic error possibility by 
comparing L = 440 h' 1 Mpc, 220 h' 1 Mpc, and 110 h' 1 Mpc 
boxes, with the power in each case averaged over 9 realiza- 
tions with different seeds. Remarkably, there is no signifi- 
cant sign of systematic error, even in the L — 110 /i -1 Mpc 
calculation (the error in the absolute power is larger). The 
maximum disagreement in Figure |3J < 2% between the 
L = 440 /i -1 Mpc and L = 220 h' 1 Mpc simulations at 
2 = 1.5, appears not to be strictly a boxsize effect at all, 
but rather a coupling between boxsize and limited particle 
density (the particle noise appears to whiten more quickly 
in the larger box). For these comparisons we have matched 



the particle density and force resolution between the two 
box sizes, i.e., we compare (440, 192, 768) to (220, 96, 384), 
and (220, 192, 768) to (110, 96, 384). 

3.6 Time steps 

For our main grid we used ~ 210 time steps to evolve the 
simulations (~ 120 for 0.0 < z < 1.5). Reducing this to 
~ 80 (~ 30) leads to the < 1% change shown in Figure [TU1 
It appears that we could accelerate our grid calculation by 
a factor of a few by relaxing our standard time step restric- 
tions. 

3.7 Power spectrum computation 

We generally use an (N = 1024) 3 grid for the power spec- 
trum computation, with a simple correction for the CIC 
smoothing, which we see in Figure 1111 is essentially exact 
(better than 0.5% as tested using measurements wit h dif- 
ferent TV) out to ~ 0.7 fcNyq, where fcNyq = tyN/L (see I Jin J 
(2004) for a discussion of aliasing). Note that if we are look- 
ing at a simple ratio of raw power in two simulations of 
the same size, the CIC smoothing correction we apply has 
no effect, because it is multiplicative (assuming we are not 
subtracting shot-noise). Figure HTI shows that a grid spacing 
Ax ~ 0.2 Mpc is sufficient to introduce essentially no 
error into our computation at k < 10 /iMpc -1 , independent 
of the particle density. 
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Figure 10. Effect of time step size on the fractional effect of 
w on the non-linear power. The black/solid (blue/short-dashed) 
line shows at z = (z = 1.5) the ratio of f w —~ o.s(^) computed 
from (220, 192, 768) simulations with ~ 80 time steps (~ 30 for 
0.0 < z < 1.5) to ~ 210 (~ 120) time steps. The red/long-dashed 
(green/dotted) line shows the same comparison for (110, 192, 768) 
simulations. 



4 SIMULATION GRID AND FITTING 
FORMULA 

Our limited quantitative goal in th is paper is to prov ide 
a module that can be grafted onto ISmith et alJ i2003l) to 
account for varying w. 

We found in §3 that to reduce finite particle density and 
force resolution errors to 2-3%, we need resolution equivalent 
to (110, 192, 768). An L = 110 h" 1 Mpc box is sufficient to 
achieve percent level systematic accuracy, and ~ 3% statis- 
tical precision, i.e., (110, 192, 768) simulations are sufficient 
to essentially solve the problem of computing the effect of 
ra to a few percent, especially if we average over a few dif- 
ferent realizations of the initial conditions, and considering 
that the maximum errors are generally near 10 /i Mpc -1 , 
where baryons probably limit our precision anyway. The 
code we describe in this section is based on a grid of (110, 
192, 768) simulations, averaged over four realizations of the 
initial conditions for each grid point. At this point it would 
be straightforward to perform a grid of larger simulations 
to meet the numerical requirements more comfortably, but 
given the generally limited scope of this paper, we have de- 
ferred this to future work. 

Our grid of models is motivated by the idea of Taylor 
expanding the dependence of -Pnl(w)/-Pnl(u> = — 1) on the 
other parameters around a central model. The central model 
and positive and negative variations are as = 0.897 ± 0.097, 
fi m = 0.281±0.070, Q b = 0.0462±0.0052, h = 0.710±0.066, 
and n = 0. 980 ± 0.065 ( motiva ted by the best fit and 3a 
errors from ISeliak et alJ i2005l) ) To be clear: we are only 
varying these parameters individually, not in combinations, 



Figure 11. Change in effect of w, fw——o.s(^) = 
P w —_Q c t (k)/P w — _i, when we change the resolution of the grid 
that we use to compute the power spectrum for L = 110 h~ 1 Mpc 
simulations. The denominator is always power computed using 
(N = 1024) 3 cells. For P = 192, M = 768, black (solid) show the 
change for N = 512 and red (dotted) N = 256, while for P = 96, 
M = 384 we use blue (long-dashed) and green (short-dashed) to 
show the same two iV's, respectively. The thick lines show 2 = 0, 
thin lines show z = 1.5. The vertical cyan (dotted) lines mark 
(128,256)tt/L. 



i.e., the grid does not have 3^ points. For the two most 
important parameters, fl m and as, we add the four possible 
joint variations to the grid. For each variation of the non-ui 
parameters, we run models with w =-0.5, -0.75, -1.0, and 
-1.5. We extract the power spectrum at z =1.5, 1.0, 0.5, 
0.25, and 0.0. We do not advocate this kind of grid for a more 
general simulation project, where a set of sim ulations guided 
by C MB constraints should be most efficient dWhite fc Vald 

It wou l d pro bably be straightforward to extend 
ISmith et alJ J20Q3>) by modifying their fitting functions 
/i(fira) to depend on w, but we take the less sophisti- 
cated approach of describing the change in power relative 
to w — — 1 by a multi-polynomial function of the cosmologi- 
cal parameters. If p is the vector of cosmological parameters, 
which we take to include z, the formula for the correction 
factor is 



In 



Pnl(u>)/-D 2 0) 



(fc,p) = 



(3) 



Pnl(w = -1)/D2(w = -1) 

(N p Ni \ 
(*) , 
i=li/ i= J 

where N p is the number of cosmological parameters, Ni 
is the order of polynomial to use for each of them, and 
A v1V2 u 3 ...vn (M) are coefficients to be determined by a least- 
squares fit to simulations. D(w) is the linear growth factor, 
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normalized to 1 at z = (we divide out the growth factor in 
Equation 0| to remove the relatively trivial linear evolution 
with redshift). Note that this formula includes many cross- 
terms for which we do not have simulations in our grid - their 
coefficients are set to zero when the fit is performed using 
singular value decomposition (we use Equation0]because it 
is easy to write down and implement in code, and extends 
automatically when additional simulations become avail- 
able). We measure Pnl(A;) from the simulations in bands 
spaced by Alog 10 k = 0.1, and determine a separate set of 
As for each band. Once the As are determined, Pnl for any 
model is computed by simply plugging the desired parame- 
ters into Equation 0| and multiplying the result by the non- 
linear power in the corresponding w = — 1 model, and the 
appropriate linear growth factors. Our code quantifying the 
results, with an example showing how to use it, can be found 
at http://www.cita.utoronto.ca/~pmcdonal/code.html un- 
der the name "wcorrector." This code is only tested for 
k < 10/iMpc -1 . The user who needs to integrate to higher 
k should absorb the uncertainty in the extrapolation into 
the uncertainty they are assuming from baryon effects. 

Note that we have not separated the dependence of the 
non-linear power on fl m through its effect on past non-linear 
growth from its effect through the transfer function (except 
in Figure I^J. Separating these would probably be useful in 
the future, so that formulas can be used with arbitrary lin- 
ear power spectrum without fear that spurious effects of Q m 
will be generated. However, because we have used fully ac - 
curate transfer functions from "lingers" feertschingerl200ll) . 
the coupling between the two influences of f2 ro can have no 
effect on calculations within the standard ACDM model. Ad- 
ditionally, we showed that the effect of changing the power 
spectrum shape on f m at fixed as is negligible, so it can 
probably be safely assumed that the effect of fi m that we 
see is coming through the past growth. 

The reader may ask whether a more thoughtful, theoret- 
ically motivated fitting formula could be more general and 
efficient. This is possible. On the other hand, a drawback 
of such formulas is that they trap the user into a limited 
functional form for the various dependences. When a care- 
fully constructed fitting formula is found not to fit to the 
required level of precision, it is not generally straightforward 
to extend it. It is simple to extend Equation [I] to fit any set 
of simula t ions. Furthermore, as we saw in the case of the 
iMa et alJ <1999l) formula, it is easy in these cases to be de- 
ceived into thinking that they are more generally applicable 
than they really are, i.e., to think of them as having gen- 
uine predictive power rather than simply being a method of 
interpolation between simulation results. Equation 0] makes 
the interpolatory nature of these fitting formulas explicit. 
Ultimately, a hybrid approach may be optimal. A physically 
motivated fitting formula could be used to remove much of 
the gross parameter dependence, with a general formula like 
Equation0]used to make corrections for the imperfections in 
the fitting formula. The hope would be that this would allow 
a sparser calibration grid than would otherwise be needed. 



5 CONCLUSIONS 

We have isolated the effect of changing w on the non-linear 
power spectrum of mass density fluctuations by compar- 



ing simulations with identical linear theory density fields 
at the observed redshift. We focused on this definition of 
the effect of w (i.e., fixed linear theory power and other 
parameters at z = 0) simply because it has not been care- 
fully considered in the past, and this comple ments predic- 
tive fo rmulas calibrated only for w = — 1 (e.g.. ISmith et ail 
(2003)). The change in power relative to w = —1 is <10% 
for k < 1/iMpc" 1 (for -0.5 < w < -1.5), but rises to 12% 
by k — 10/iMpc" 1 in a model with w — —0.75, and ~ 33% 
for w = —0.5 (at z = 0). Among the other cosmological 
parameters, the size of the effect is primarily sensitive to 
the dark energy fraction, i.e., Q m in flat models. The power 
spectrum normalization (i.e., as) also has a small effect, but 
the slope/shape of the power spectrum are irrelevant (as 
represented by varying n, fit, and h at fixed as). 

Figure con firms the accuracy of the formula of 
ISmith et all lf2003l) for the dependence on fl m in ACDM 
models, at least once it has been modified to use high accu- 
racy transfer functions. 

We provide a simple code 

(http:/ /www. cita.utoronto.ca/~pmcdonal/code. html ) 
quantifying the effect of w as a function of k, z, w, Q m , 
as, n, h, and Qb, to be used as a correction to Pnl(A;) 
calculations accurate at w = —1. The dependence of the 
power spectrum on w should be accurate to a few percent 
for fc<10 /iMpc -1 . Our quantitative results may be useful 
for making more realistic projections of the future potential 
to measure w by methods sensitive to the non-linear power 
(primarily weak lensing). Our code will be appropriate 
for forecasts of parameter measurements from future large 
data sets, or parameter determinations using data sets that 
at least include WMAP. It should be used cautiously for 
fits to limited amounts of weak lensing data alone, since 
we do not cover extreme parameter values (however, by 
construction the code will only be unreliable in regions of 
parameter space that are ruled out by other observations). 

Our ambitions have been quite limited in this paper. 
We have not tried to determine the absolute power in the 
central model because that is much more difficult to simu- 
late precisely than the fractional changes we studied here, 
requiring both larger box sizes and higher resolution. Much 
of the error caused by limitations in the simulations cancels 
when we take ratios of power spectra, a fact that should 
make future construction of high precision grids of simula- 
tion predictions easier. 
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